* ADL DGP 

// Basics
clear all
set matsize 11000

parallel setclusters 8

capture program drop adl
program define adl, rclass

// Monte Carlo Basics
version 12.1
syntax [, obs(integer 1) rho_e(real 0.5) rho_x(real 0.5) ]
drop _all
*args obs rho_e rho_x
set obs `obs'

// Time Series Basics
egen time = fill(0 1 2)
tsset time

// Data Generating Process

* Unit root error
gen e = rnormal()

* AR(1) predictor
gen x = rnormal()
replace x = `rho_x'*L.x + rnormal() if time > 0

* Generate outcome
scalar beta = 2
scalar gamma1 = 1.5
gen y = rnormal()
replace y = `rho_e'*L.y + beta*x + gamma1*L.x + e if time > 0

// Model Estimation
* ADL (1,1)
reg y L.y x L.x
return scalar phi_adl = _b[L1.y]
scalar phi_adl_test_t1 = (_b[L1.y] - `rho_e') / _se[L1.y]
return scalar phi_adl_t1 = 2*ttail(e(df_r),abs(phi_adl_test_t1)) < 0.05
scalar phi_adl_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi_adl_power = 2*ttail(e(df_r),abs(phi_adl_test_power)) < 0.05
return scalar b1_adl = _b[x]
scalar b1_adl_test_t1 = (_b[x] - beta) / _se[x]
return scalar b1_adl_t1 = 2*ttail(e(df_r),abs(b1_adl_test_t1)) < 0.05
scalar b1_adl_test_power = (_b[x] - 0) / _se[x]
return scalar b1_adl_power = 2*ttail(e(df_r),abs(b1_adl_test_power)) < 0.05
return scalar b2_adl = _b[L1.x]
scalar b2_adl_test_t1 = (_b[L1.x] - gamma1) / _se[L1.x]
return scalar b2_adl_t1 = 2*ttail(e(df_r),abs(b2_adl_test_t1)) < 0.05
scalar b2_adl_test_power = (_b[L1.x] - 0) / _se[L1.x]
return scalar b2_adl_power = 2*ttail(e(df_r),abs(b2_adl_test_power)) < 0.05
return scalar lrm_adl = (_b[x] + _b[L1.x]) /(1-_b[L1.y])
nlcom ((_b[x] + _b[L1.x]) /(1-_b[L1.y])) * (1e+3)
matrix define var_lrm_adl = r(V)
scalar var_lrm_adl_se  = sqrt(var_lrm_adl[1,1]) / (1e+3)
scalar lrm_adl_test_t1 = (((_b[x] + _b[L1.x]) /(1-_b[L1.y])) - ((beta + gamma1)/(1-`rho_e'))) /var_lrm_adl_se
return scalar lrm_adl_t1 = 2*ttail(e(df_r),abs(lrm_adl_test_t1)) < 0.05
scalar lrm_adl_test_power = (((_b[x] + _b[L1.x]) /(1-_b[L1.y])) - 0) /var_lrm_adl_se
return scalar lrm_adl_power = 2*ttail(e(df_r),abs(lrm_adl_test_power)) < 0.05

* GECM
reg D.y L.y D.x L.x
return scalar phi_gecm = _b[L1.y]
local gecm_true_phi = `rho_e' - 1
scalar phi_gecm_test_t1 = (_b[L1.y] - `gecm_true_phi') / _se[L1.y]
return scalar phi_gecm_t1 = 2*ttail(e(df_r),abs(phi_gecm_test_t1)) < 0.05
scalar phi_gecm_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi_gecm_power = 2*ttail(e(df_r),abs(phi_gecm_test_power)) < 0.05
return scalar b1_gecm  = _b[D1.x]
scalar b1_gecm_test_t1= (_b[D1.x] - beta) / _se[D1.x]
return scalar b1_gecm_t1 = 2*ttail(e(df_r),abs(b1_gecm_test_t1)) < 0.05
scalar b1_gecm_test_power = (_b[D1.x] - 0) / _se[D1.x]
return scalar b1_gecm_power = 2*ttail(e(df_r),abs(b1_gecm_test_power)) < 0.05
return scalar b2_gecm  = _b[L1.x]
scalar b2_gecm_test_t1 = (_b[L1.x] - (beta + gamma1)) / _se[L1.x]
return scalar b2_gecm_t1 = 2*ttail(e(df_r),abs(b2_gecm_test_t1)) < 0.05
scalar b2_gecm_test_power = (_b[L1.x] - 0) / _se[L1.x]
return scalar b2_gecm_power = 2*ttail(e(df_r),abs(b2_gecm_test_power)) < 0.05
return scalar lrm_gecm = -_b[L1.x] / _b[L1.y]
nlcom (-_b[L1.x] / _b[L1.y])* (1e+3)
matrix define var_lrm_gecm = r(V)
scalar var_lrm_gecm_se  = sqrt(var_lrm_gecm[1,1]) / (1e+3)
scalar lrm_gecm_test_t1 = ((-_b[L1.x] / _b[L1.y]) - ((beta + gamma1)/(1-`rho_e'))) /var_lrm_gecm_se
return scalar lrm_gecm_t1 = 2*ttail(e(df_r),abs(lrm_gecm_test_t1)) < 0.05
scalar lrm_gecm_test_power = ((-_b[L1.x] / _b[L1.y]) - 0) /var_lrm_gecm_se
return scalar lrm_gecm_power = 2*ttail(e(df_r),abs(lrm_gecm_test_power)) < 0.05

* ADL (2,2)
reg y L.y L2.y x L.x L2.x

return scalar phi1_adl22 = _b[L1.y]
scalar phi1_adl22_test_t1 = (_b[L1.y] - `rho_e') / _se[L1.y]
return scalar phi1_adl22_t1 = 2*ttail(e(df_r),abs(phi1_adl22_test_t1)) < 0.05
scalar phi1_adl22_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi1_adl22_power = 2*ttail(e(df_r),abs(phi1_adl22_test_power)) < 0.05
return scalar b1_adl22 = _b[x]
return scalar phi2_adl22 = _b[L2.y]
scalar phi2_adl22_test = (_b[L2.y] - 0) / _se[L2.y]
return scalar phi2_adl22_t1 = 2*ttail(e(df_r),abs(phi2_adl22_test)) < 0.05
scalar b1_adl22_test_t1 = (_b[x] - beta) / _se[x]
return scalar b1_adl22_t1 = 2*ttail(e(df_r),abs(b1_adl22_test_t1)) < 0.05
scalar b1_adl22_test_power = (_b[x] - 0) / _se[x]
return scalar b1_adl22_power = 2*ttail(e(df_r),abs(b1_adl22_test_power)) < 0.05
return scalar b2_adl22 = _b[L1.x]
scalar b2_adl22_test_t1 = (_b[L1.x] - gamma1 ) / _se[L1.x]
return scalar b2_adl22_t1 = 2*ttail(e(df_r),abs(b2_adl22_test_t1)) < 0.05
scalar b2_adl22_test_power = (_b[L1.x] - 0) / _se[L1.x]
return scalar b2_adl22_power = 2*ttail(e(df_r),abs(b2_adl22_test_power)) < 0.05
return scalar b3_adl22 = _b[L2.x]
scalar b3_adl22_test = (_b[L2.x] - 0) / _se[L2.x]
return scalar b3_adl22_t1 = 2*ttail(e(df_r),abs(b3_adl22_test)) < 0.05
return scalar lrm_adl22 = (_b[x] + _b[L1.x] + _b[L2.x]) /(1-_b[L1.y]-_b[L2.y])
nlcom ((_b[x] + _b[L1.x] + _b[L2.x]) /(1-_b[L1.y]-_b[L2.y])) * (1e+3)
matrix define var_lrm_adl22 = r(V)
scalar var_lrm_adl22_se  = sqrt(var_lrm_adl22[1,1]) / (1e+3)
scalar lrm_adl22_test_t1 = (((_b[x] + _b[L1.x] + _b[L2.x]) /(1-_b[L1.y]-_b[L2.y])) - ((beta + gamma1)/(1-`rho_e'))) /var_lrm_adl22_se
return scalar lrm_adl22_t1 = 2*ttail(e(df_r),abs(lrm_adl22_test_t1)) < 0.05
scalar lrm_adl22_test_power = (((_b[x] + _b[L1.x] + _b[L2.x]) /(1-_b[L1.y]-_b[L2.y])) - 0) /var_lrm_adl22_se
return scalar lrm_adl22_power = 2*ttail(e(df_r),abs(lrm_adl22_test_power)) < 0.05


* ADL (3,3)
reg y L.y L2.y L3.y x L.x L2.x L3.x

return scalar phi1_adl33 = _b[L1.y]
scalar phi1_adl33_test_t1 = (_b[L1.y] - `rho_e') / _se[L1.y]
return scalar phi1_adl33_t1 = 2*ttail(e(df_r),abs(phi1_adl33_test_t1)) < 0.05
scalar phi1_adl33_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi1_adl33_power = 2*ttail(e(df_r),abs(phi1_adl33_test_power)) < 0.05
return scalar phi2_adl33 = _b[L2.y]
scalar phi2_adl33_test = (_b[L2.y] - 0) / _se[L2.y]
return scalar phi2_adl33_t1 = 2*ttail(e(df_r),abs(phi2_adl33_test)) < 0.05
return scalar phi3_adl33 = _b[L3.y]
scalar phi3_adl33_test = (_b[L3.y] - 0) / _se[L3.y]
return scalar phi3_adl33_t1 = 2*ttail(e(df_r),abs(phi3_adl33_test)) < 0.05
return scalar b1_adl33 = _b[x]
scalar b1_adl33_test_t1 = (_b[x] - beta) / _se[x]
return scalar b1_adl33_t1 = 2*ttail(e(df_r),abs(b1_adl33_test_t1)) < 0.05
scalar b1_adl33_test_power = (_b[x] - 0) / _se[x]
return scalar b1_adl33_power = 2*ttail(e(df_r),abs(b1_adl33_test_power)) < 0.05
return scalar b2_adl33 = _b[L1.x]
scalar b2_adl33_test_t1 = (_b[L1.x] - gamma1 ) / _se[L1.x]
return scalar b2_adl33_t1 = 2*ttail(e(df_r),abs(b2_adl33_test_t1)) < 0.05
scalar b2_adl33_test_power = (_b[L1.x] - 0) / _se[L1.x]
return scalar b2_adl33_power = 2*ttail(e(df_r),abs(b2_adl33_test_power)) < 0.05
return scalar b3_adl33 = _b[L2.x]
scalar b3_adl33_test = (_b[L2.x] - 0) / _se[L2.x]
return scalar b3_adl33_t1 = 2*ttail(e(df_r),abs(b3_adl33_test)) < 0.05
return scalar b4_adl33 = _b[L3.x]
scalar b4_adl33_test = (_b[L3.x] - 0) / _se[L3.x]
return scalar b4_adl33_t1 = 2*ttail(e(df_r),abs(b4_adl33_test)) < 0.05
return scalar lrm_adl33 = (_b[x] + _b[L1.x] + _b[L2.x] + _b[L3.x]) /(1-_b[L1.y]-_b[L2.y]-_b[L3.y])
nlcom ((_b[x] + _b[L1.x]+ _b[L2.x] + _b[L3.x]) /(1-_b[L1.y]-_b[L2.y]-_b[L3.y])) * (1e+3)
matrix define var_lrm_adl33 = r(V)
scalar var_lrm_adl33_se  = sqrt(var_lrm_adl33[1,1]) / (1e+3)
scalar lrm_adl33_test_t1 = (((_b[x] + _b[L1.x]+ _b[L2.x] + _b[L3.x]) /(1-_b[L1.y]-_b[L2.y]-_b[L3.y])) - ((beta + gamma1))/(1-`rho_e')) /var_lrm_adl33_se
return scalar lrm_adl33_t1 = 2*ttail(e(df_r),abs(lrm_adl33_test_t1)) < 0.05
scalar lrm_adl33_test_power = (((_b[x] + _b[L1.x]+ _b[L2.x] + _b[L3.x]) /(1-_b[L1.y]-_b[L2.y]-_b[L3.y])) - 0) /var_lrm_adl33_se
return scalar lrm_adl33_power = 2*ttail(e(df_r),abs(lrm_adl33_test_power)) < 0.05

* F-test, ADL(2,2)
reg y L.y L2.y x L.x L2.x
test L.y + L2.y = `rho_e'
return scalar ftest_y_adl22_t1 = r(p) < 0.05
reg y L.y L2.y x L.x L2.x
test L.y + L2.y = 0
return scalar ftest_y_adl22_power = r(p) < 0.05
reg y L.y L2.y x L.x L2.x
test x + L.x + L2.x = beta + gamma1
return scalar ftest_x_adl22_t1 = r(p) < 0.05
reg y L.y L2.y x L.x L2.x
test x + L.x + L2.x = 0
return scalar ftest_x_adl22_power = r(p) < 0.05

* F-test, ADL(3,3)
reg y L.y L2.y L3.y x L.x L2.x L3.x
test L.y + L2.y + L3.y = `rho_e'
return scalar ftest_y_adl33_t1 = r(p) < 0.05
reg y L.y L2.y L3.y x L.x L2.x L3.x
test L.y + L2.y + L3.y = 0
return scalar ftest_y_adl33_power = r(p) < 0.05
reg y L.y L2.y L3.y x L.x L2.x L3.x
test x + L.x + L2.x + L3.x = beta + gamma1
return scalar ftest_x_adl33_t1 = r(p) < 0.05
reg y L.y L2.y L3.y x L.x L2.x L3.x
test x + L.x + L2.x + L3.x = 0
return scalar ftest_x_adl33_power = r(p) < 0.05

end
